Análisis exploratorio de los datos#
Contextualización#
import pandas as pd
import matplotlib.pyplot as plt
import plotly.express as px
import seaborn as sns
import numpy as np
import statsmodels.tsa.api as smtsa
from statsmodels.tsa.stattools import adfuller
from statsmodels.graphics.tsaplots import plot_acf
from statsmodels.graphics.tsaplots import plot_pacf
from statsmodels.tsa.seasonal import seasonal_decompose
from multiprocessing import Pool, cpu_count
Se conoce como PM10 a aquellas partículas sólidas o líquidas de diferente composición que se encuentran dispersas en la atmósfera. Estas pueden ser causadas por consecuencias humanas o naturales. Las características más importantes de estas partículas tienen un diámetro aerodinámico menor que 10 \(\mu m\). Al ser estas partículas tan pequeñas, estas pueden ser inhaladas por nuestro sistema respiratorio, por esto se les conoce a estas partículas como fracción respirable.
Es de interés el estudio de la presencia de estas partículas, puesto que generan efectos adversos sobre la salud, como asma agravada, función pulmonar reducida, irritación en las vías respiratorias, hasta muerte prematura en personas con enfermedades cardíacas o pulmonares.
En este proyecto serán usados datos que han sido recolectados en las estaciones de muestreo de la calidad del aire por el DAGMA (Departamento Administrativo de Gestión del Medio Ambiente) de la ciudad de Cali.
Estaciones de monitoreo.#
La Flora
ERA - Obrero
Transitoria EDB-Navarro
Base Aérea
Pance
Univalle
Compartir
La Ermita
Cañaveralejo
Análisis#
data = pd.read_csv("C:/Users/fonta/OneDrive/Escritorio/data_2017_2022.csv", sep=";")
data.tail()
| fecha | estacion | variable | tipo | medicion | |
|---|---|---|---|---|---|
| 4084056 | 2022-12-31T23:59:59Z | era_obrero | velocidad_viento | val | 0,4000000000000000 |
| 4084057 | 2022-12-31T23:59:59Z | canaveralejo | velocidad_viento | val | NaN |
| 4084058 | 2022-12-31T23:59:59Z | canaveralejo | direccion_viento | val | NaN |
| 4084059 | 2022-12-31T23:59:59Z | canaveralejo | lluvia | val | NaN |
| 4084060 | 2022-12-31T23:59:59Z | canaveralejo | radiacion_solar | val | NaN |
data.shape
(4084061, 5)
Nuestro dataset cuenta con 4.084.061 de registros inicialmente, pero luego de unos procesos estos registros se reducirán.Estos registros fueron tomados desde diferentes bases de medición y observados en el tiempo.
data.fecha.unique()
array(['2017-01-01T00:59:59Z', '2017-01-01T01:59:59Z',
'2017-01-01T02:59:59Z', ..., '2022-12-31T21:59:59Z',
'2022-12-31T22:59:59Z', '2022-12-31T23:59:59Z'], dtype=object)
Como primer análisis de nuestro conjunto de datos, vemos que se tienen registros desde el año 2017 hasta 2022, estos registros estan tomados por cada hora de cada día durante estos 6 años.
Variables únicas
data.variable.unique()
array(['pm25', 'so2', 'lluvia', 'o3', 'temperatura', 'pm10', 'humedad',
'velocidad_viento', 'direccion_viento', 'radiacion_solar',
'presion', 'no2', 'uv-pm', 'black_carbon', 'h2s',
'temperatura_10_m'], dtype=object)
Nuestro dataset tiene un total de 16 variables, pero recordemos que el interés de este proyecto es predecir la comtamicación por la particula del Pm10, así que mas adelante trabajaremos solo con los registros correspondientes de esta variable.
Este dataset viene organizado de una manera no tan adecuada para su estudio, es por eso que es necesario hacer un pivotado, con el fin de poder tener acceso mejor a los datos de interés.
Pivotado del conjunto de datos#
data1=data.pivot(index=['fecha', 'estacion'],columns='variable',values='medicion')
data1= data1.reset_index()
data1.head(5)
| variable | fecha | estacion | black_carbon | direccion_viento | h2s | humedad | lluvia | no2 | o3 | pm10 | pm25 | presion | radiacion_solar | so2 | temperatura | temperatura_10_m | uv-pm | velocidad_viento |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 2017-01-01T00:59:59Z | base_aerea | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN |
| 1 | 2017-01-01T00:59:59Z | canaveralejo | NaN | NaN | NaN | NaN | NaN | NaN | NaN | 51,7000000000000028 | NaN | NaN | NaN | 0,3142645629907114 | NaN | NaN | NaN | NaN |
| 2 | 2017-01-01T00:59:59Z | compartir | NaN | 208,5999999999999943 | NaN | 75,9000000000000057 | 0,0000000000000000 | NaN | 8,4562520617694883 | 221,0000000000000000 | 161,0000000000000000 | 680,5000000000000000 | 0,0000000000000000 | NaN | 25,3000000000000007 | NaN | NaN | 1,3999999999999999 |
| 3 | 2017-01-01T00:59:59Z | era_obrero | NaN | NaN | NaN | NaN | 0,0000000000000000 | NaN | 16,1473212223579807 | 195,8000000000000114 | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN |
| 4 | 2017-01-01T00:59:59Z | ermita | NaN | NaN | NaN | NaN | NaN | NaN | NaN | 61,8999999999999986 | NaN | NaN | NaN | 2,1474745137698608 | NaN | NaN | NaN | NaN |
data1.shape
(473256, 18)
Luego de los procesos anteriores vemos que ahora tenemos un nuevo dataset, esta vez con las variables organizadas adecuadamente y con su respectiva medición, así obtenemos un conjunto final de 473.256 registros y 18 variables.
Veamos las proporciones de datos faltantes.
missing_data = data1.isna().sum()
print(missing_data)
variable
fecha 0
estacion 0
black_carbon 464636
direccion_viento 268837
h2s 457903
humedad 242369
lluvia 212585
no2 473255
o3 271320
pm10 238688
pm25 298710
presion 264075
radiacion_solar 251972
so2 341076
temperatura 206050
temperatura_10_m 459510
uv-pm 464636
velocidad_viento 267851
dtype: int64
Según lo anterior vemos cuantos valores faltantes por cada variable, en este caso es de nuestro interés solo la variable PM10, que en esta caso cuenta con 238.688 registros faltantes, mas adelante veremos como trataremos estos datos faltantes.
Es importante ajustar el formato de algunas de las variables por ejemplo la fecha:
data1['fecha'] = pd.to_datetime(data1['fecha'], utc=True)
data1['fecha']
0 2017-01-01 00:59:59+00:00
1 2017-01-01 00:59:59+00:00
2 2017-01-01 00:59:59+00:00
3 2017-01-01 00:59:59+00:00
4 2017-01-01 00:59:59+00:00
...
473251 2022-12-31 23:59:59+00:00
473252 2022-12-31 23:59:59+00:00
473253 2022-12-31 23:59:59+00:00
473254 2022-12-31 23:59:59+00:00
473255 2022-12-31 23:59:59+00:00
Name: fecha, Length: 473256, dtype: datetime64[ns, UTC]
A continuación debemos ajustar los caracteres de puntuación en la variable de interés.
data1['pm10'] = data1['pm10'].str.replace(',', '.', regex=False)
data1['pm10']= pd.to_numeric(data1['pm10'], errors='coerce')
data1['pm10']
0 NaN
1 51.7
2 221.0
3 195.8
4 61.9
...
473251 NaN
473252 NaN
473253 NaN
473254 NaN
473255 NaN
Name: pm10, Length: 473256, dtype: float64
Puesto que la cantidad de datos faltantes es muy grande y con el fin de hacer un estudio de manera específica, se harán unos filtros por cada estación y se trabajará para este proyecto la estación que tenga menos datos faltantes.
Filtro por estación#
baseaerea = data1[(data1["estacion"]=="base_aerea") ]
compartir = data1[(data1["estacion"]=="compartir") ]
univalle = data1[(data1["estacion"]=="univalle") ]
canaveralejo = data1[(data1["estacion"]=="canaveralejo") ]
ermita = data1[(data1["estacion"]=="ermita") ]
eraobrero = data1[(data1["estacion"]=="era_obrero") ]
flora = data1[(data1["estacion"]=="flora") ]
pance = data1[(data1["estacion"]=="pance") ]
transitoria = data1[(data1["estacion"]=="transitoria") ]
baseaerea_na = baseaerea["pm10"].isna().sum()
compartir_na = compartir["pm10"].isna().sum()
univalle_na = univalle["pm10"].isna().sum()
canaveralejo_na = canaveralejo["pm10"].isna().sum()
ermita_na = ermita["pm10"].isna().sum()
eraobrero_na = eraobrero["pm10"].isna().sum()
flora_na = flora["pm10"].isna().sum()
pance_na = pance["pm10"].isna().sum()
transitoria_na = transitoria["pm10"].isna().sum()
values_pm10=[baseaerea_na,compartir_na,univalle_na,canaveralejo_na,ermita_na,eraobrero_na,flora_na,pance_na,transitoria_na]
values_pm10=sorted(values_pm10)
print(values_pm10)
print('')
print("NA en pm10 para base_aerea:", baseaerea_na)
print("NA en pm10 para base_aerea:", compartir_na)
print("NA en pm10 para univalle:", univalle_na)
print("NA en pm10 para canaveralejo:", canaveralejo_na)
print("NA en pm10 para ermita:", ermita_na)
print("NA en pm10 para era_obrero:", eraobrero_na)
print("NA en pm10 para flora:", flora_na)
print("NA en pm10 para base_aerea:", pance_na)
print("NA en pm10 para transitoria:", transitoria_na)
[13410, 14301, 16819, 16986, 17459, 21497, 40761, 44871, 52584]
NA en pm10 para base_aerea: 52584
NA en pm10 para base_aerea: 14301
NA en pm10 para univalle: 40761
NA en pm10 para canaveralejo: 21497
NA en pm10 para ermita: 13410
NA en pm10 para era_obrero: 17459
NA en pm10 para flora: 16986
NA en pm10 para base_aerea: 16819
NA en pm10 para transitoria: 44871
Al hacer los respectivos filtros y cuentas vemos que la estación con menos valores faltantes es la estación de ERMITA así que de ahora en adelante se trabajará solo con esta base de datos.
Ermita#
ermita.info()
<class 'pandas.core.frame.DataFrame'>
Index: 52584 entries, 4 to 473251
Data columns (total 18 columns):
# Column Non-Null Count Dtype
--- ------ -------------- -----
0 fecha 52584 non-null datetime64[ns, UTC]
1 estacion 52584 non-null object
2 black_carbon 0 non-null object
3 direccion_viento 0 non-null object
4 h2s 0 non-null object
5 humedad 26997 non-null object
6 lluvia 0 non-null object
7 no2 0 non-null object
8 o3 0 non-null object
9 pm10 39174 non-null float64
10 pm25 13878 non-null object
11 presion 22874 non-null object
12 radiacion_solar 0 non-null object
13 so2 33476 non-null object
14 temperatura 32687 non-null object
15 temperatura_10_m 0 non-null object
16 uv-pm 0 non-null object
17 velocidad_viento 0 non-null object
dtypes: datetime64[ns, UTC](1), float64(1), object(16)
memory usage: 7.6+ MB
len(ermita)
52584
A partir de lo anterior tenemos que la base de datos de la estación de Ermita cuenta con 52.584 registros de los cuales la variable PM10 tiene 39.174 no nulos.
Se creará un dataset con solo nuestras variables de interés que son la fecha y la medición del PM10.
Serie de tiempo de PM10 en la estación ermita#
serie_pm10= ermita[['fecha', 'pm10']].reset_index()
serie_pm10.drop(columns=['index'],inplace=True)
serie_pm10
| variable | fecha | pm10 |
|---|---|---|
| 0 | 2017-01-01 00:59:59+00:00 | 61.9 |
| 1 | 2017-01-01 01:59:59+00:00 | 101.0 |
| 2 | 2017-01-01 02:59:59+00:00 | 75.6 |
| 3 | 2017-01-01 03:59:59+00:00 | 56.3 |
| 4 | 2017-01-01 04:59:59+00:00 | 107.1 |
| ... | ... | ... |
| 52579 | 2022-12-31 19:59:59+00:00 | NaN |
| 52580 | 2022-12-31 20:59:59+00:00 | NaN |
| 52581 | 2022-12-31 21:59:59+00:00 | NaN |
| 52582 | 2022-12-31 22:59:59+00:00 | NaN |
| 52583 | 2022-12-31 23:59:59+00:00 | NaN |
52584 rows × 2 columns
Veamos la distribucion de la variable:
sns.histplot(serie_pm10['pm10'],kde=True, color='#473C8B')
plt.title('Distribución de la variable PM10')
plt.xlabel('')
plt.ylabel('')
plt.show()
Vemos que la distribución de la variable es asimétrica a la derecha, aun así presenta la forma de una distribución normal.
Intentemos graficar la serie de tiempo para la estación ermita.
fig = px.line(serie_pm10,x='fecha', y='pm10', title='Serie de tiempo PM10 2017-2022')
fig.update_traces(line_color='#473C8B')
fig.show()
Al ver la gráfica de la serie de tiempo no se diferencia mucho lo que es la tendencia o el patrón de la serie de tiempo, pero lo que si podemos distinguir claramente son esos registros nulos, podemos ver que la mayor cantidad de datos faltantes se encuentran en el año 2021, con esta imagen podemos ir pensando en como se hará la imputación de estos datos faltantes.
Veamos la proporción de datos faltantes en la serie de tiempo.
print("Los valores faltantes del Pm10 en la estación ermita son en total: ",ermita_na)
print("Proporción de los valores faltantes: ", ermita_na/len(ermita))
Los valores faltantes del Pm10 en la estación ermita son en total: 13410
Proporción de los valores faltantes: 0.25502053856686446
Como muestra la gráfica de la serie y como comentamos anteriormente, la mayoria de los datos faltantes estan a partir de 2021, por esta razón se hará el análisis e imputación de los datos faltantes dividiendo el dataset. La metodología que se usará es la siguiente: Se evaluara e implementará un modelo con los datos desde 2017 hasta 2019 y se imputarán los datos faltantes con un modelo ARIMA, cuando se cumpla la imputación y se verifique que fue correcta, se tomaran todos estos datos como entrenamiento para poder hacer la imputación de los datos faltantes que estan a partir de 2020 hasta 2022. Veamos la implementación.
Serie de tiempo 2017-2019#
A continuación se crea la serie de tiempo del pm10 en el periodo de 2017 a 2019.
ermita17_19=serie_pm10[serie_pm10['fecha'].dt.year <=2019]
ermita17_19
| variable | fecha | pm10 |
|---|---|---|
| 0 | 2017-01-01 00:59:59+00:00 | 61.9 |
| 1 | 2017-01-01 01:59:59+00:00 | 101.0 |
| 2 | 2017-01-01 02:59:59+00:00 | 75.6 |
| 3 | 2017-01-01 03:59:59+00:00 | 56.3 |
| 4 | 2017-01-01 04:59:59+00:00 | 107.1 |
| ... | ... | ... |
| 26275 | 2019-12-31 19:59:59+00:00 | 38.3 |
| 26276 | 2019-12-31 20:59:59+00:00 | 37.3 |
| 26277 | 2019-12-31 21:59:59+00:00 | 40.0 |
| 26278 | 2019-12-31 22:59:59+00:00 | 28.9 |
| 26279 | 2019-12-31 23:59:59+00:00 | 35.0 |
26280 rows × 2 columns
Veamos como es el comportamiento de la serie con los siguientes análisis.
Comportamiento de la serie#
Inicialmente omitimos los valores NaN.
serie_pm10_wna= ermita17_19.dropna(subset=['pm10']).reset_index(drop=True)
serie_pm10_wna
| variable | fecha | pm10 |
|---|---|---|
| 0 | 2017-01-01 00:59:59+00:00 | 61.9 |
| 1 | 2017-01-01 01:59:59+00:00 | 101.0 |
| 2 | 2017-01-01 02:59:59+00:00 | 75.6 |
| 3 | 2017-01-01 03:59:59+00:00 | 56.3 |
| 4 | 2017-01-01 04:59:59+00:00 | 107.1 |
| ... | ... | ... |
| 24023 | 2019-12-31 19:59:59+00:00 | 38.3 |
| 24024 | 2019-12-31 20:59:59+00:00 | 37.3 |
| 24025 | 2019-12-31 21:59:59+00:00 | 40.0 |
| 24026 | 2019-12-31 22:59:59+00:00 | 28.9 |
| 24027 | 2019-12-31 23:59:59+00:00 | 35.0 |
24028 rows × 2 columns
Como bien es sabido es importante en el estudio de las series de tiempo verificar como es su gráfico de autocorrelación, el cual muestra la correlación que tiene un registros con sus registros en el pasado, veamos.
plt.rcParams.update({'figure.figsize': (10, 4)})
plot_acf(serie_pm10_wna["pm10"], lags=100)
plt.tight_layout()
plt.show()
Si analisamos el gráfico de autocorrelación vemos que existe un patrón entre las correlaciones de los valores anteriores de la serie, esto es un indicio para pensar en que estamos trabajando con una serie estacionaria o incluso cíclica, vemos que los dos primeros valores de autocorrelación son cernas a 1, es decir que hay una dependencia secuencial en los datos.
Es importante notar que a medida que aumentan los lags la correlacion va disminuyendo.
Así mismo como se estudia la autocorrelación, tambien es importante estudiar el gráfico de autocorrelación parcial, el cual muestra la correlación que tiene un registros con lo observado antes en el pasado eliminando la influencia de los registros intermedios.
plt.rcParams.update({'figure.figsize': (10, 4)})
plot_pacf(serie_pm10_wna["pm10"], lags=100)
plt.tight_layout()
plt.show()
Por parte del gráfico de autocorrelación parcial vemos un cambio abrupto entre las primeras correlaciones, indicando así que cada observación esta fuertemente influenciada por la observacioón inmediata anterior a diferencia de las observaciones mas distantes en el pasado.
Se realizará la prueba de Duckey-Fuller aumentada para probar de forma analítica que la serie es estacionaria.
result = adfuller(serie_pm10_wna["pm10"])
print('ADF Statistic: %f' % result[0])
print('p-value: %f' % result[1])
ADF Statistic: -15.342534
p-value: 0.000000
Si recordamos la hipotesis nula de la prueba de Duckey-Fuller, esta se plantea como $\(H_0= \textit{La serie no es estacionaria}\)$ Con un p-valore = 0 se rechaza la hipotesis nula para cualquier nivel de significacia, por tanto se prueba que la serie es estacionaria.
Puesto que se confirma que la serie de tiempo es estacionaria, se implementarán uno de los modelos clásicos de las series de tiempo para hacer la imputación de los datos faltantes en este caso se hará uso de un modelo ARIMA o mejor dicho un modelo ARMA puesto que la serie ya es estacionaria, a continuación se presentan la hiperparametrización del modelo y mas adelante la predicción de los datos faltantes.
Modelo ARMA para la imputación de datos faltantes.#
Se hará una busqueda de los mejores parametros para nuestro modelo arma:
aicVal=[]
for ari in range(1, 3):
for maj in range(1,3):
arma_obj = smtsa.ARIMA(serie_pm10_wna["pm10"].tolist(), order=(ari, 0, maj)).fit()
aicVal.append([ari, maj, arma_obj.aic])
dfAIC = pd.DataFrame(aicVal, columns=['AR(p)', 'MA(q)', 'AIC'])
dfAIC
| AR(p) | MA(q) | AIC | |
|---|---|---|---|
| 0 | 1 | 1 | 199288.621318 |
| 1 | 1 | 2 | 199275.645535 |
| 2 | 2 | 1 | 199258.876717 |
| 3 | 2 | 2 | 199259.913104 |
Luego de hacer lo que se podría llamar un GridSearch para los parametros del modelo ARMA, se creó un data set donde se registran los scores de cada uno de los modelos con las diferentes combinaciones de parametros. Este proceso arrojó los siguientes resultados.
dfAIC.nsmallest(n=1, columns="AIC")
| AR(p) | MA(q) | AIC | |
|---|---|---|---|
| 2 | 2 | 1 | 199258.876717 |
Vemos que al hacer la busqueda de los mejores parametros el proceso arroja como mejores parametros \(p=2\) y \(q=1\) es decir que el modelo que mejor se ajusta a los datos de la base de Ermita es ARMA(2,1).
Implementación de la imputación:
Primeramente necesitamos el df que se querrá imputar, es decir la serie con los valores falatantes.
to_imp=ermita17_19.copy()
to_imp
| variable | fecha | pm10 |
|---|---|---|
| 0 | 2017-01-01 00:59:59+00:00 | 61.9 |
| 1 | 2017-01-01 01:59:59+00:00 | 101.0 |
| 2 | 2017-01-01 02:59:59+00:00 | 75.6 |
| 3 | 2017-01-01 03:59:59+00:00 | 56.3 |
| 4 | 2017-01-01 04:59:59+00:00 | 107.1 |
| ... | ... | ... |
| 26275 | 2019-12-31 19:59:59+00:00 | 38.3 |
| 26276 | 2019-12-31 20:59:59+00:00 | 37.3 |
| 26277 | 2019-12-31 21:59:59+00:00 | 40.0 |
| 26278 | 2019-12-31 22:59:59+00:00 | 28.9 |
| 26279 | 2019-12-31 23:59:59+00:00 | 35.0 |
26280 rows × 2 columns
Cantidad de valores a imputar.
to_imp["pm10"].isna().sum()
2252
Vemos que se busca imputar 2.252 observaciones.
Por otro lado se necesita la serie sin los NaN para poder entrenar el modelo ARMA y posteriormente usar las predicciones del modelo para imputar los datos faltantes.
entrenamiento=serie_pm10_wna.copy()
entrenamiento
| variable | fecha | pm10 |
|---|---|---|
| 0 | 2017-01-01 00:59:59+00:00 | 61.9 |
| 1 | 2017-01-01 01:59:59+00:00 | 101.0 |
| 2 | 2017-01-01 02:59:59+00:00 | 75.6 |
| 3 | 2017-01-01 03:59:59+00:00 | 56.3 |
| 4 | 2017-01-01 04:59:59+00:00 | 107.1 |
| ... | ... | ... |
| 24023 | 2019-12-31 19:59:59+00:00 | 38.3 |
| 24024 | 2019-12-31 20:59:59+00:00 | 37.3 |
| 24025 | 2019-12-31 21:59:59+00:00 | 40.0 |
| 24026 | 2019-12-31 22:59:59+00:00 | 28.9 |
| 24027 | 2019-12-31 23:59:59+00:00 | 35.0 |
24028 rows × 2 columns
Vemos que se usarán 24.028 observaciones para entrenar el modelo que nos ayudará a predecir los datos faltantes.
Se necesita de igual forma los indices donde se encuentran los valore faltantes con el fin de rellenarlos con las predicciones del modelo.
missing_indices = to_imp[to_imp['pm10'].isnull()].index
print(len(missing_indices))
missing_indices
2252
Index([ 168, 169, 170, 730, 731, 1021, 1046, 1047, 1048, 1049,
...
24224, 24742, 24947, 24948, 24949, 25427, 25428, 25429, 25760, 25761],
dtype='int64', length=2252)
Ajuste del modelo#
En el paso siguiente se hará el entrenamiento del modelo con los parametros escogidos anteriormente, luego se haran las predicciónes y se extraerán los valores en los indices donde se encuentran los valores faltantes.
arma_obj_finp = smtsa.ARIMA(entrenamiento["pm10"].tolist(), order=(2, 0, 1)).fit()
predicciones = arma_obj_finp.predict(start=0, end=len(to_imp) - 1)
for idx in missing_indices:
to_imp.at[idx, 'pm10'] = predicciones[idx]
Rectifiquemos que se haya hecho la imputación de forma correcta.
to_imp["pm10"].isna().sum()
0
Recordemos como se veía la serie con los valores faltantes a continuación.
fig = px.line(ermita17_19,x='fecha', y='pm10', title='Serie de tiempo PM10 2017-2022')
fig.update_traces(line_color='#473C8B')
fig.show()
Y ahora veamos como se ve la serie con los valores imputados.
fig = px.line(to_imp,x='fecha', y='pm10', title='Serie de tiempo PM10 2017-2022')
fig.update_traces(line_color='#473C8B')
fig.show()
Con el gráfico anterior vemos que la imputación se hizo de forma correcta aparentemente, revisemos los supuestos del modelo para confirmar de manera analítica como se comportó el modelo.
Revisión de supuestos#
arma_obj_finp.plot_diagnostics(figsize=(10,7))
plt.show()
Los residuos, que representan las diferencias entre los valores observados y los predichos por el modelo, no presentan patrones claros de autocorrelación y su varianza parece ser constante a lo largo del tiempo. Estos hallazgos sugieren que los residuos son esencialmente ruido aleatorio, ademas si notamos el histrograma y el QQ plot, parece ser que los residuos se distribuyen normales y no guardan una correlación entre ellos.
Ahora veamos que pasa en la serie de tiempo a partir del 2020. La idea es usar los datos que anteriormente imputamos para poder imputar los siguientes, puesto que en esta parte del dataframe se encuentra la mayor parte de los datos faltantes.
#Este conjunto de datos será ahora nuestro conjunto de entrenamiento para el modelo.
entrenamiento_2=to_imp.copy()
entrenamiento_2
| variable | fecha | pm10 |
|---|---|---|
| 0 | 2017-01-01 00:59:59+00:00 | 61.9 |
| 1 | 2017-01-01 01:59:59+00:00 | 101.0 |
| 2 | 2017-01-01 02:59:59+00:00 | 75.6 |
| 3 | 2017-01-01 03:59:59+00:00 | 56.3 |
| 4 | 2017-01-01 04:59:59+00:00 | 107.1 |
| ... | ... | ... |
| 26275 | 2019-12-31 19:59:59+00:00 | 38.3 |
| 26276 | 2019-12-31 20:59:59+00:00 | 37.3 |
| 26277 | 2019-12-31 21:59:59+00:00 | 40.0 |
| 26278 | 2019-12-31 22:59:59+00:00 | 28.9 |
| 26279 | 2019-12-31 23:59:59+00:00 | 35.0 |
26280 rows × 2 columns
Serie de tiempo 2020-2022#
ermita20_22=serie_pm10[serie_pm10['fecha'].dt.year >=2020].reset_index(drop=True)
ermita20_22
| variable | fecha | pm10 |
|---|---|---|
| 0 | 2020-01-01 00:59:59+00:00 | NaN |
| 1 | 2020-01-01 01:59:59+00:00 | NaN |
| 2 | 2020-01-01 02:59:59+00:00 | NaN |
| 3 | 2020-01-01 03:59:59+00:00 | NaN |
| 4 | 2020-01-01 04:59:59+00:00 | NaN |
| ... | ... | ... |
| 26299 | 2022-12-31 19:59:59+00:00 | NaN |
| 26300 | 2022-12-31 20:59:59+00:00 | NaN |
| 26301 | 2022-12-31 21:59:59+00:00 | NaN |
| 26302 | 2022-12-31 22:59:59+00:00 | NaN |
| 26303 | 2022-12-31 23:59:59+00:00 | NaN |
26304 rows × 2 columns
Recordemos al inicio del estudio de la serie de tiempo que mencionamos que a partir del 2020 esta la mayor cantidad de datos faltantes.
Comportamiento de la serie#
Iniciamos omitiendo los NaN para el estudio de los aspectos importantes de la serie.
serie_pm10_wna= ermita20_22.dropna(subset=['pm10']).reset_index(drop=True)
serie_pm10_wna
| variable | fecha | pm10 |
|---|---|---|
| 0 | 2020-01-01 10:59:59+00:00 | 13.8 |
| 1 | 2020-01-01 11:59:59+00:00 | 9.4 |
| 2 | 2020-01-01 12:59:59+00:00 | 11.5 |
| 3 | 2020-01-01 13:59:59+00:00 | 15.1 |
| 4 | 2020-01-01 14:59:59+00:00 | 19.6 |
| ... | ... | ... |
| 15141 | 2022-11-01 05:59:59+00:00 | 57.7 |
| 15142 | 2022-11-01 06:59:59+00:00 | 77.5 |
| 15143 | 2022-11-01 07:59:59+00:00 | 77.2 |
| 15144 | 2022-11-01 08:59:59+00:00 | 89.4 |
| 15145 | 2022-11-01 23:59:59+00:00 | 27.7 |
15146 rows × 2 columns
Vemos que el total de los valores en la serie desde 2020 que no son faltantes 15.146 es decir que mas de la mitad de los datos son faltantes.
Pasemos a revisar el gráfico de autocorrelación.
plt.rcParams.update({'figure.figsize': (10,4)})
plot_acf(serie_pm10_wna["pm10"], lags=100)
plt.tight_layout()
plt.show()
Vemos que el gráfico de autocorrelacion de esta parte de la serie es muy similar, muestra de igual forma un patron de correlaciones, lo que nos lleva a pensar a que de igual forma esta serie es estacionaria, es importante tambien notar que a medida que nos dezplazamos en el pasado las correlaciones van disminuyendo.
plt.rcParams.update({'figure.figsize': (8,4)})
plot_pacf(serie_pm10_wna["pm10"], lags=50)
plt.tight_layout()
plt.show()
Y como era de esperarse en la gráfica de la autocorrelacion parcial se ve el mismo patrón que obtuvimos en la primera parte de la serie, vemos que de igual manera las dos primeras correlaciones indican que esta serie esta muy correlacionada con sus dos observaciones inmediatas del pasado.
De igual forma para confirmar, se realiza a continuación la prueba se Duckey-Fuller.
result = adfuller(serie_pm10_wna["pm10"])
print('ADF Statistic: %f' % result[0])
print('p-value: %f' % result[1])
ADF Statistic: -12.381327
p-value: 0.000000
Mostrando así lo anterior, que la serie de tiempo desde 2020 hasta 2022 es estacionaria tambien ya que se obtuvo un p-valor=0 rechazando la hipotesis nula bajo todos los niveles de significancia.
Según lo anterior tambien es correcto implentar un modelo ARMA, ya que la serie es estacionaria y este modelo es el indicado para llevar a cabo la imputación.
Modelo ARMA para la imputación de datos faltantes#
Hiperparametrización de los parametros p y q para el modelo. Recordemos que utilizaremos los valores que imputamos en al seccion anterior del dataset, es por eso que debemos hacer una concatenación entre los valores desde 2017 hasta 2019 y con los valores no nulos de esta sección.
serie_pm10_wna=pd.concat([entrenamiento_2, serie_pm10_wna])
len(serie_pm10_wna)
41426
aicVal=[]
for ari in range(1, 4):
for maj in range(1,4):
arma_obj = smtsa.ARIMA(serie_pm10_wna["pm10"].tolist(), order=(ari, 0, maj)).fit()
aicVal.append([ari, maj, arma_obj.aic])
C:\Users\fonta\miniconda3\envs\ml_venv\lib\site-packages\statsmodels\tsa\statespace\sarimax.py:966: UserWarning:
Non-stationary starting autoregressive parameters found. Using zeros as starting parameters.
C:\Users\fonta\miniconda3\envs\ml_venv\lib\site-packages\statsmodels\tsa\statespace\sarimax.py:978: UserWarning:
Non-invertible starting MA parameters found. Using zeros as starting parameters.
dfAIC = pd.DataFrame(aicVal, columns=['AR(p)', 'MA(q)', 'AIC'])
dfAIC
| AR(p) | MA(q) | AIC | |
|---|---|---|---|
| 0 | 1 | 1 | 342142.877387 |
| 1 | 1 | 2 | 342124.701464 |
| 2 | 1 | 3 | 342100.221717 |
| 3 | 2 | 1 | 342111.514901 |
| 4 | 2 | 2 | 342111.050623 |
| 5 | 2 | 3 | 342100.977151 |
| 6 | 3 | 1 | 342111.242819 |
| 7 | 3 | 2 | 342112.585103 |
| 8 | 3 | 3 | 342102.073159 |
dfAIC.nsmallest(n=1, columns="AIC")
| AR(p) | MA(q) | AIC | |
|---|---|---|---|
| 2 | 1 | 3 | 342100.221717 |
Este proceso nos arroja que los mejores parametros para ajustar nuestro modelo ARMA es \(p=1\) y \(q=3\).
Implementación de la imputación:
Se tiene el dataset con los valores faltantes que se quieren imputar.
to_impk=ermita20_22.copy()
to_impk
| variable | fecha | pm10 |
|---|---|---|
| 0 | 2020-01-01 00:59:59+00:00 | NaN |
| 1 | 2020-01-01 01:59:59+00:00 | NaN |
| 2 | 2020-01-01 02:59:59+00:00 | NaN |
| 3 | 2020-01-01 03:59:59+00:00 | NaN |
| 4 | 2020-01-01 04:59:59+00:00 | NaN |
| ... | ... | ... |
| 26299 | 2022-12-31 19:59:59+00:00 | NaN |
| 26300 | 2022-12-31 20:59:59+00:00 | NaN |
| 26301 | 2022-12-31 21:59:59+00:00 | NaN |
| 26302 | 2022-12-31 22:59:59+00:00 | NaN |
| 26303 | 2022-12-31 23:59:59+00:00 | NaN |
26304 rows × 2 columns
Total de valores que se busca imputar:
to_impk["pm10"].isna().sum()
11158
Vemos que el total de valores que se quiere imputar son 11.158 datos.
Recordemos cual será nuestro conjunto de entrenamiento.
serie_pm10_wna
| variable | fecha | pm10 |
|---|---|---|
| 0 | 2017-01-01 00:59:59+00:00 | 61.9 |
| 1 | 2017-01-01 01:59:59+00:00 | 101.0 |
| 2 | 2017-01-01 02:59:59+00:00 | 75.6 |
| 3 | 2017-01-01 03:59:59+00:00 | 56.3 |
| 4 | 2017-01-01 04:59:59+00:00 | 107.1 |
| ... | ... | ... |
| 15141 | 2022-11-01 05:59:59+00:00 | 57.7 |
| 15142 | 2022-11-01 06:59:59+00:00 | 77.5 |
| 15143 | 2022-11-01 07:59:59+00:00 | 77.2 |
| 15144 | 2022-11-01 08:59:59+00:00 | 89.4 |
| 15145 | 2022-11-01 23:59:59+00:00 | 27.7 |
41426 rows × 2 columns
El siguiente paso es gurdar donde se encuentran los valores faltantes, con el fin de poder hacer la imputación de manera correcta.
missing_indices = to_impk[to_impk['pm10'].isnull()].index
print(len(missing_indices))
missing_indices
11158
Index([ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9,
...
26294, 26295, 26296, 26297, 26298, 26299, 26300, 26301, 26302, 26303],
dtype='int64', length=11158)
Ajustando el modelo.#
arma_obj_finp = smtsa.ARIMA(serie_pm10_wna["pm10"].tolist(), order=(1, 0, 3)).fit()
predicciones = arma_obj_finp.predict(start=0, end=len(to_impk) - 1)
for idx in missing_indices:
to_impk.at[idx, 'pm10'] = predicciones[idx]
Revisamos que la imputación se haya hecho de la forma correcta.
to_impk["pm10"].isna().sum()
0
fig = px.line(ermita20_22,x='fecha', y='pm10', title='Serie de tiempo PM10 2017-2022')
fig.update_traces(line_color='#473C8B')
fig.show()
Veamos la gráfica de la serie de tiempo desde el 2020 hasta 2022 con la presencia de valores faltantes.
Y ahora vemos la gráfica con los valores imputados.
fig = px.line(to_impk,x='fecha', y='pm10', title='Serie de tiempo PM10 2017-2022')
fig.update_traces(line_color='#473C8B')
fig.show()
Vemos que aparentemente el modelo logró imputar de manera correcta los valores faltantes de esta segunda sección de la serie de tiempo, revisemos los supuestos del modelo para corroborar esta afirmación.
Revisión de supuestos#
arma_obj_finp.plot_diagnostics(figsize=(12,7))
plt.show()
Si verificamos los supuestos del modelo podemos ver que los residuos se distribuyen normales y que no tiene una correlación entre ellos lo cual es importante, puesto que indica que los errores que existen son aleatorios. Por tanto la imputación realizada por el modelo ARIMA (1,0,3).
Serie final#
serie_pm10_imputada=pd.concat([to_imp,to_impk])
serie_pm10_imputada.to_csv('seriepm10.csv',index=False)
fig = px.line(serie_pm10_imputada,x='fecha', y='pm10', title='Serie de tiempo PM10 2017-2022')
fig.update_traces(line_color='#473C8B')
fig.show()
Vemos de forma gráfica como fue el resultado de la imputación, veamos si luego de esta se mantuvo la ditribución de la variable:
fig, axs = plt.subplots(1, 2, figsize=(12, 5))
# Histograma 1
sns.histplot(serie_pm10['pm10'], kde=True, color='#473C8B', ax=axs[0])
axs[0].set_title('Distribución original')
axs[0].set_xlabel('')
axs[0].set_ylabel('')
# Histograma 2
sns.histplot(serie_pm10_imputada['pm10'], kde=True, color='#473C8B', ax=axs[1])
axs[1].set_title('Distribución luego de la imputación')
axs[1].set_xlabel('')
axs[1].set_ylabel('')
# Mostrar el gráfico
plt.tight_layout()
plt.show()
Vemos que la imputación no alteró la distribución de la variable que es lo que se esperaba.
Analisemos los gráficos de autocorrelación y autocorrelación parcial de la serie final.
plt.rcParams.update({'figure.figsize': (10,4)})
plot_acf(serie_pm10_imputada["pm10"], lags=100)
plt.tight_layout()
plt.show()
plt.rcParams.update({'figure.figsize': (8,4)})
plot_pacf(serie_pm10_imputada["pm10"], lags=50)
plt.tight_layout()
plt.show()
Con esto podemos darnos un indicio de que la serie de tiempo del PM10 con los datos imputados, tambien es estacionaria., lo cual es lo epserado, puesto que la serie original tambien lo era, es decir, no se alteró el patrón de la serie de tiempo.
Con este análisis exploratorio de datos podemos pasar a implementar los diferentes modelos, con el fin de decidir cual es el que mejor funciona para nuestro conjunto de datos. Pasemos a la modelación.